Wind trajectory
library(GeoPressureR)
library(tidyverse)
library(leaflet)
library(leaflet.extras)
library(raster)
library(dplyr)
library(ggplot2)
library(plotly)
knitr::opts_chunk$set(echo = FALSE)
load(paste0("../data/1_pressure//", params$gdl_id, "_pressure_prob.Rdata"))
load(paste0("../data/2_light/", params$gdl_id, "_light_prob.Rdata"))
load(paste0("../data/3_static/", params$gdl_id, "_static_prob.Rdata"))
# load(paste0("../data/4_basic_graph/", params$gdl_id, "_basic_graph.Rdata"))
load(paste0("../data/5_wind_graph/", params$gdl_id, "_wind_graph.Rdata"))
load(paste0("../data/5_wind_graph/", params$gdl_id, "_grl.Rdata"))
col <- rep(RColorBrewer::brewer.pal(8, "Dark2"), times = ceiling(max(tag$sta$sta_id) / 8))
Altitudes are computed based on pressure measurement of the geolocation, corrected based on the assumed location of the shortest path. This correction accounts therefore for the natural variation of pressure as estimated by ERA-5. The vertical lines indicate the sunrise (dashed) and sunset (solid).
p <- ggplot() +
# geom_line(data = tag$pressure, aes(x = date, y = obs), colour = "grey") +
geom_line(data = do.call("rbind", shortest_path_timeserie), aes(x = date, y = altitude)) +
geom_line(data = do.call("rbind", shortest_path_timeserie) %>% filter(sta_id > 0), aes(x = date, y = altitude, col = factor(sta_id))) +
geom_vline(data = twl, aes(xintercept = twilight, linetype = ifelse(rise, "dashed", "solid"), color="grey"), lwd=0.1) +
theme_bw() +
scale_colour_manual(values = col) +
scale_y_continuous(name = "Altitude (m)")
ggplotly(p, dynamicTicks = T) %>% layout(showlegend = F)
file <- paste0("figure_print/wintering_location/wintering_location_",params$gdl_id,".png")
if(file.exists(file)){
knitr::include_graphics(file)
}
tmp <- lapply(pressure_prob, function(x) {
mt <- metadata(x)
df <- data.frame(
start = mt$temporal_extent[1],
end = mt$temporal_extent[2],
sta_id = mt$sta_id
)
})
tmp2 <- do.call("rbind", tmp)
sim_lat <- as.data.frame(t(path_sim$lat)) %>%
mutate(sta_id = path_sim$sta_id) %>%
pivot_longer(-c(sta_id)) %>%
left_join(tmp2,by="sta_id")
sim_lat_p <- sim_lat %>%
filter(sta_id==max(sta_id)) %>%
mutate(start=end) %>%
rbind(sim_lat)
sp_lat <- as.data.frame(shortest_path) %>% left_join(tmp2,by="sta_id")
sp_lat_p <- sp_lat %>%
filter(sta_id==max(sta_id)) %>%
mutate(start=end) %>%
rbind(sp_lat)
p <- ggplot() +
geom_step(data=sim_lat_p, aes(x=start, y=value, group=name), alpha=.07) +
geom_point(data=sp_lat_p, aes(x=start, y=lat)) +
xlab('Date') +
ylab('Latitude') +
theme_light()
ggplotly(p, dynamicTicks = T)
The large circles indicates the shortest path (overall most likely trajectory) estimated by the graph approach. The size is proportional to the duration of stay. The small dots and grey lines represents 10 possible trajeectories of the bird according to the model.
Click on the full-screen mode button on the top-left of the map to see more details on the map.
sta_duration <- unlist(lapply(static_prob_marginal, function(x) {
as.numeric(difftime(metadata(x)$temporal_extent[2], metadata(x)$temporal_extent[1], units = "days"))
}))
pal <- colorFactor(col, as.factor(seq_len(length(col))))
m <- leaflet(width = "100%") %>%
addProviderTiles(providers$Stamen.TerrainBackground) %>%
addFullscreenControl() %>%
addPolylines(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = "#808080", weight = 3) %>%
addCircles(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = pal(factor(shortest_path$sta_id, levels = tag$sta$sta_id)), weight = sta_duration^(0.3) * 10)
for (i in seq_len(nrow(path_sim$lon))) {
m <- m %>%
addPolylines(lng = path_sim$lon[i, ], lat = path_sim$lat[i, ], opacity = 0.5, weight = 1, color = "#808080") %>%
addCircles(lng = path_sim$lon[i, ], lat = path_sim$lat[i, ], opacity = .7, weight = 1, color = pal(factor(shortest_path$sta_id, levels = tag$sta$sta_id)))
}
m
The marginal probability map estimate the overall probability of position at each stationary period regardless of the trajectory taken by the bird. It is the most useful quantification of the uncertainty of the position of the bird.
li_s <- list()
l <- leaflet(width = "100%") %>%
addProviderTiles(providers$Stamen.TerrainBackground) %>%
addFullscreenControl()
for (i_r in seq_len(length(static_prob_marginal))) {
i_s <- metadata(static_prob_marginal[[i_r]])$sta_id
info <- metadata(static_prob_marginal[[i_r]])$temporal_extent
info_str <- paste0(i_s, " | ", info[1], "->", info[2])
li_s <- append(li_s, info_str)
l <- l %>%
addRasterImage(static_prob_marginal[[i_r]], colors = "OrRd", opacity = 0.8, group = info_str) %>%
addCircles(lng = shortest_path$lon[i_s], lat = shortest_path$lat[i_s], opacity = 1, color = "#000", weight = 10, group = info_str)
}
l %>%
addLayersControl(
overlayGroups = li_s,
options = layersControlOptions(collapsed = FALSE)
) %>%
hideGroup(tail(li_s, length(li_s) - 1))
fun_marker_color <- function(norm){
if (norm < 20){
"darkpurple"
} else if (norm < 35){
"darkblue"
} else if (norm < 50){
"lightblue"
} else if (norm < 60){
"lightgreen"
} else if (norm < 80){
"yellow"
} else if (norm < 100){
"lightred"
} else {
"darkred"
}
}
fun_NSEW <- function(angle){
angle <- angle %% (pi* 2)
angle <- angle*180/pi
if (angle < 45/2){
"E"
} else if (angle < 45*3/2){
"NE"
} else if (angle < 45*5/2){
"N"
} else if (angle < 45*7/2){
"NW"
} else if (angle < 45*9/2){
"W"
} else if (angle < 45*11/2){
"SW"
} else if (angle < 45*13/2){
"S"
}else if (angle < 45*15/2){
"SE"
} else {
"E"
}
}
sta_duration <- unlist(lapply(static_prob_marginal,function(x){as.numeric(difftime(metadata(x)$temporal_extent[2],metadata(x)$temporal_extent[1],units="days"))}))
m <-leaflet(width = "100%") %>%
addProviderTiles(providers$Stamen.TerrainBackground) %>% addFullscreenControl() %>%
addPolylines(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = "#808080", weight = 3) %>%
addCircles(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = "#000", weight = sta_duration^(0.3)*10)
for (i_s in seq_len(grl$sz[3]-1)){
if (grl$flight_duration[i_s]>5){
edge <- which(grl$s == shortest_path$id[i_s] & grl$t == shortest_path$id[i_s+1])
label = paste0( i_s,': ', grl$flight[[i_s]]$start, " - ", grl$flight[[i_s]]$end, "<br>",
"F. dur.: ", round(grl$flight_duration[i_s]), ' h <br>',
"GS: ", round(abs(grl$gs[edge])), ' km/h, ',fun_NSEW(Arg(grl$gs[edge])),'<br>',
"WS: ", round(abs(grl$ws[edge])), ' km/h, ',fun_NSEW(Arg(grl$ws[edge])),'<br>',
"AS: ", round(abs(grl$as[edge])), ' km/h, ',fun_NSEW(Arg(grl$as[edge])),'<br>')
iconArrow <- makeAwesomeIcon(icon = "arrow-up",
library = "fa",
iconColor = "#FFF",
iconRotate = (90 - Arg(grl$ws[edge])/pi*180) %% 360,
squareMarker = TRUE,
markerColor = fun_marker_color(abs(grl$ws[edge])))
m <- m %>% addAwesomeMarkers(lng = (shortest_path$lon[i_s] + shortest_path$lon[i_s+1])/2,
lat = (shortest_path$lat[i_s] + shortest_path$lat[i_s+1])/2,
icon = iconArrow, popup = label)
}
}
m
edge <- t(graph_path2edge(path_sim$id, grl))
nj <- ncol(edge)
nsta <- ncol(path_sim$lon)
speed_df <- data.frame(
as = abs(grl$as[edge]),
gs = abs(grl$gs[edge]),
ws = abs(grl$ws[edge]),
sta_id_s = rep(head(grl$sta_id,-1), nj),
sta_id_t = rep(tail(grl$sta_id,-1), nj),
flight_duration = rep(head(grl$flight_duration,-1), nj),
dist = geosphere::distGeo(
cbind(as.vector(t(path_sim$lon[,1:nsta-1])), as.vector(t(path_sim$lat[,1:nsta-1]))),
cbind(as.vector(t(path_sim$lon[,2:nsta])), as.vector(t(path_sim$lat[,2:nsta])))
) / 1000
) %>% mutate(
name = paste(sta_id_s,sta_id_t, sep="-")
)
plot1 <- ggplot(speed_df, aes(reorder(name, sta_id_s), gs)) + geom_boxplot() + theme_bw() +scale_x_discrete(name = "")
plot2 <- ggplot(speed_df, aes(reorder(name, sta_id_s), ws)) + geom_boxplot() + theme_bw() +scale_x_discrete(name = "")
plot3 <- ggplot(speed_df, aes(reorder(name, sta_id_s), as)) + geom_boxplot() + theme_bw() +scale_x_discrete(name = "")
plot4 <- ggplot(speed_df, aes(reorder(name, sta_id_s), flight_duration)) + geom_point() + theme_bw() +scale_x_discrete(name = "")
subplot(ggplotly(plot1), ggplotly(plot2), ggplotly(plot3), ggplotly(plot4), nrows=4, titleY=T)
alt_df = do.call("rbind", shortest_path_timeserie) %>%
arrange(date) %>%
mutate(
sta_id_s = cummax(sta_id),
sta_id_t = sta_id_s+1
) %>%
filter(sta_id == 0 & sta_id_s > 0 ) %>%
group_by(sta_id_s, sta_id_t) %>%
summarise(
alt_min = min(altitude),
alt_max = max(altitude),
alt_mean = mean(altitude),
alt_med = median(altitude),
)
trans_df <- speed_df %>%
group_by(sta_id_s,sta_id_t,flight_duration) %>%
summarise(
as_m = mean(as),
as_s = sd(as),
gs_m = mean(gs),
gs_s = sd(gs),
ws_m = mean(ws),
ws_s = sd(ws),
dist_m = mean(dist),
dist_s = sd(dist)
) %>%
left_join(alt_df)
trans_df %>% kable()
| sta_id_s | sta_id_t | flight_duration | as_m | as_s | gs_m | gs_s | ws_m | ws_s | dist_m | dist_s | alt_min | alt_max | alt_mean | alt_med |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 2 | 3.4166667 | 34.57011 | 17.494468 | 47.68620 | 21.711983 | 18.003929 | 3.8248656 | 162.92785 | 74.18261 | 157.25331 | 1060.0105 | 654.1848 | 604.5139 |
| 2 | 3 | 7.3333333 | 37.18249 | 15.983395 | 44.39281 | 16.646386 | 11.408376 | 5.5214162 | 325.54730 | 122.07350 | 240.07618 | 985.7225 | 595.2773 | 547.5521 |
| 3 | 4 | 7.9166667 | 40.58991 | 15.977018 | 56.06959 | 18.303820 | 19.826997 | 5.5771796 | 443.88426 | 144.90524 | -230.29992 | 1404.8480 | 298.0325 | 214.8520 |
| 4 | 5 | 4.9166667 | 33.81206 | 17.370380 | 34.80925 | 16.913821 | 4.123158 | 2.5650512 | 171.14546 | 83.15962 | 22.06228 | 1241.4214 | 487.8901 | 506.4481 |
| 5 | 6 | 9.8333333 | 41.63452 | 11.463145 | 50.63115 | 16.145250 | 11.780202 | 2.6895206 | 497.87295 | 158.76162 | 9.47344 | 1670.3913 | 679.8558 | 587.9694 |
| 6 | 7 | 9.4166667 | 36.63912 | 13.033330 | 36.98587 | 12.525579 | 12.530809 | 5.0572737 | 348.28359 | 117.94921 | 127.30825 | 3458.0496 | 1838.6165 | 1684.6414 |
| 7 | 9 | 28.4166667 | 55.26767 | 17.586456 | 64.81036 | 14.110518 | 11.531872 | 3.7032344 | 1841.69430 | 400.97388 | NA | NA | NA | NA |
| 9 | 10 | 8.5000000 | 62.58095 | 19.842185 | 71.09937 | 18.510293 | 16.076736 | 2.4089642 | 604.34461 | 157.33749 | 935.15268 | 2205.3068 | 1804.8507 | 1854.0280 |
| 10 | 11 | 8.9166667 | 41.04600 | 17.261752 | 48.03535 | 18.698191 | 20.129450 | 4.1477964 | 428.31520 | 166.72553 | 658.94857 | 3100.9877 | 1691.4406 | 1707.8348 |
| 11 | 12 | 5.2500000 | 41.59498 | 21.239332 | 40.42210 | 19.901930 | 3.833050 | 1.4844871 | 212.21603 | 104.48513 | 295.45738 | 1527.9584 | 1134.8679 | 1336.0644 |
| 12 | 13 | 7.1666667 | 44.52462 | 17.710289 | 38.32131 | 17.694810 | 8.304144 | 1.8732288 | 274.63602 | 126.81281 | 359.33122 | 1215.7742 | 706.7753 | 650.7422 |
| 13 | 14 | 0.8333333 | 35.96519 | 16.702519 | 36.62206 | 20.024884 | 9.981134 | 3.0560772 | 30.51839 | 16.68740 | 136.05294 | 219.7510 | 177.9019 | 177.9019 |
| 14 | 16 | 12.6666667 | 12.42522 | 6.712675 | 12.08427 | 7.529162 | 3.545367 | 0.9657721 | 153.06746 | 95.36939 | NA | NA | NA | NA |
| 16 | 17 | 3.7500000 | 32.95460 | 11.883280 | 32.75673 | 11.478415 | 2.103072 | 1.1820068 | 122.83775 | 43.04406 | 220.08353 | 411.3244 | 317.3691 | 321.5518 |
| 17 | 18 | 1.5000000 | 33.64262 | 10.438046 | 39.00593 | 17.020175 | 17.288656 | 1.6683386 | 58.50889 | 25.53026 | 141.33298 | 171.3591 | 154.2831 | 150.1572 |
| 18 | 19 | 10.6666667 | 37.47178 | 11.796753 | 43.69214 | 10.826953 | 18.739869 | 2.9062068 | 466.04952 | 115.48749 | 243.59053 | 1672.0438 | 896.6158 | 895.9195 |
| 19 | 20 | 4.9166667 | 33.69398 | 14.300953 | 42.22981 | 17.339492 | 19.528778 | 4.9322794 | 207.62988 | 85.25250 | 312.52187 | 1782.4563 | 906.9879 | 825.7417 |
| 20 | 21 | 4.2500000 | 36.99632 | 14.036391 | 28.68552 | 18.504806 | 25.468923 | 4.2639659 | 121.91347 | 78.64543 | 287.69319 | 1124.9316 | 655.0057 | 669.5111 |
| 21 | 22 | 10.6666667 | 36.38271 | 17.266830 | 30.44620 | 12.582570 | 16.071009 | 7.0997736 | 324.75946 | 134.21408 | 432.44483 | 1905.1577 | 794.3203 | 628.2020 |
| 22 | 23 | 10.5000000 | 39.79834 | 16.936986 | 38.94978 | 16.247836 | 12.807723 | 5.7257082 | 408.97265 | 170.60228 | 789.05124 | 3154.8776 | 2129.4470 | 1842.2562 |
| 23 | 24 | 10.2500000 | 39.65864 | 13.748037 | 78.32582 | 15.589634 | 43.528736 | 5.2998547 | 802.83962 | 159.79375 | 1035.30185 | 3786.9226 | 3129.5720 | 3261.6373 |
| 24 | 26 | 23.5000000 | 51.53299 | 8.173692 | 113.85681 | 3.927738 | 66.987229 | 6.5230024 | 2675.63511 | 92.30185 | NA | NA | NA | NA |
| 26 | 27 | 7.9166667 | 26.75544 | 7.335790 | 31.91815 | 8.167588 | 8.769154 | 1.2363085 | 252.68532 | 64.66007 | 340.90962 | 1229.2793 | 706.1691 | 737.0251 |
| 27 | 28 | 5.7500000 | 39.82505 | 16.393176 | 57.14521 | 21.375869 | 20.570207 | 6.4380061 | 328.58494 | 122.91125 | 474.38752 | 2804.7508 | 1480.3557 | 1389.4821 |
| 28 | 29 | 6.5000000 | 61.13500 | 15.352408 | 62.17359 | 17.253713 | 25.230900 | 5.6896699 | 404.12837 | 112.14913 | 207.11052 | 2572.7330 | 964.6875 | 450.2296 |